#pragma once

#include "../../Utility/Util/utility.hpp"
#include "../../Utility/Log/Log.hpp"
#include "../../Utility/Static/Mathf.h"
#include "../../Utility/Static/Matrixf.h"
#include "../../Utility/Data/Position/CRDFactory.hpp"
#include "../../Utility/Data/Time/TimeFactory.hpp"
#include "../../Utility/Data/NavMessageData/GNSSMesData.hpp"

namespace ns_Calculate
{
    using namespace ns_Time;
    using namespace ns_Position;
    using namespace ns_GNSSMesData;
    using namespace ns_Util;
    using namespace ns_Log;

    const std::vector<std::vector<double>> Average =
        {
            {1.2769934e-3, 1.2683230e-3, 1.2465397e-3, 1.2196049e-3, 1.2045996e-3},
            {2.9153695e-3, 2.9152299e-3, 2.9288445e-3, 2.9022565e-3, 2.9024912e-3},
            {62.610505e-3, 62.837393e-3, 63.721774e-3, 63.824265e-3, 64.285455e-3}};

    const std::vector<std::vector<double>> Amplitude =
        {
            {0, 1.2709626e-5, 2.6523662e-5, 3.4000452e-5, 4.1202191e-5},
            {0, 2.1414979e-5, 3.0160779e-5, 7.2562722e-5, 11.723375e-5},
            {0, 9.0128400e-5, 4.3497037e-5, 84.795348e-5, 170.37206e-5}};

    const std::vector<double> HeightCorrection =
        {
            2.53e-5,
            5.49e-3,
            1.14e-3};

    const std::vector<std::vector<double>> WetCorr =
        {
            {5.8021897e-4, 5.6794847e-4, 5.8118019e-4, 5.9727542e-4, 6.1641693e-4},
            {1.4275268e-3, 1.5136825e-3, 1.4572752e-3, 1.5007428e-3, 1.7599082e-3},
            {4.3472961e-2, 4.6729510e-2, 4.3908931e-2, 4.4629982e-2, 5.4736032e-2}};

    class TropoCorrect
    {
    private:
        using DryWet = std::vector<double>;

        enum LatitudeAngle
        {
            Angle15,
            Angle30,
            Angle45,
            Angle60,
            Angle75
        };

        enum Coefficient
        {
            a,
            b,
            c
        };

    public:
        TropoCorrect()
        {
        }

        static double GetHopfield(const CRDFactory::CRDPtr &rr, const CRDFactory::CRDPtr &rs, const TimeFactory::TimePtr &curTime)
        {
            CRDFactory::CRDPtr _star = CRDFactory::MakeCRDPtr(CRDFactory::ValueType::TopoCentricPolarCRDSys, CartesianCRD());
            CRDFactory::SwitchTo(CRDFactory::ValueType::TopoCentricPolarCRDSys, rs, &_star, rr);
            auto star = std::dynamic_pointer_cast<TopoCentricPolarCRD>(_star);

            TimeFactory::TimePtr _curTimeOfCommonTime = TimeFactory::MakeTimePtr(TimeFactory::ValueType::CommonTimeSys);
            TimeFactory::SwitchTo(TimeFactory::ValueType::GPSTimeSys, curTime, &_curTimeOfCommonTime);
            auto curTimeOfCommonTime = std::dynamic_pointer_cast<CommonTime>(_curTimeOfCommonTime);

            CRDFactory::CRDPtr _reciver = CRDFactory::MakeCRDPtr(CRDFactory::ValueType::GeodeticCRDSys);
            CRDFactory::SwitchTo(CRDFactory::ValueType::GeodeticCRDSys, rr, &_reciver, nullptr);
            auto reciver = std::dynamic_pointer_cast<GeodeticCRD>(_reciver);

            CRDFactory::CRDPtr _Creciver = CRDFactory::MakeCRDPtr(CRDFactory::ValueType::CartesianCRDSys);
            CRDFactory::SwitchTo(CRDFactory::ValueType::CartesianCRDSys, rr, &_reciver, nullptr);
            auto Creciver = std::dynamic_pointer_cast<CartesianCRD>(_Creciver);

            double h = reciver->GetH();

            return ImproveHopfield(*star, *Creciver, h);
        }

    private:
        // 标准hopfile模型
        static double StandardHopfield(const TopoCentricPolarCRD &input, double h, int doy)
        {
            // 标准气象元素法
            double T0 = 20;
            double P0 = 1013.25;
            double RH0 = 0.5;

            double T = T0 - 0.0065 * h;
            double P = P0 * Mathf::Pow((1 - 0.0000266 * h), 5.225);
            double RH = RH0 * Mathf::Exp(-0.0006369 * h);
            double e = RH * Mathf::Exp(-37.2465 + 0.213166 * T - 0.000256908 * T * T);

            double hd = 40136 + 148.72 * (T - 273.16);
            double hw = 11000;

            double Kd = 155.2e-7 * P / T * (hd - h);
            double Kw = 155.2e-7 * 4810 / (T * T) * e * (hw - h);

            CRDFactory::CRDPtr Rr = CRDFactory::MakeCRDPtr(CRDFactory::ValueType::GeodeticCRDSys);
            CRDFactory::SwitchTo(CRDFactory::ValueType::GeodeticCRDSys,
                                 CRDFactory::MakeCRDPtr(CRDFactory::ValueType::CartesianCRDSys, input.GetCentral()), &Rr, nullptr);

            auto geoRr = std::dynamic_pointer_cast<GeodeticCRD>(Rr);

            DryWet dw = NMF(input, geoRr->GetB(), doy);

            double ans = Kd / dw[0] + Kw / dw[1];
            return ans;
        }

        // 改进Hopfile模型
        static double ImproveHopfield(const TopoCentricPolarCRD &input, const CartesianCRD &_rr, double h)
        {
            // 标准气象元素法
            double T0 = 20;
            double P0 = 1013.25;
            double RH0 = 0.5;

            double T = T0 - 0.0065 * h;
            double P = P0 * Mathf::Pow((1 - 0.0000266 * h), 5.225);
            double RH = RH0 * Mathf::Exp(-0.0006369 * h);

            // 改进hopfile模型
            double Nd = 0.776E-4;
            double Nw = 0.373;
            double r0 = Mathf::Sqrt(_rr.GetX() * _rr.GetX() + _rr.GetY() * _rr.GetY() + _rr.GetZ() * _rr.GetZ());

            double hd = 40136 + 148.72 * (T - 273.16);
            double hw = 11000;

            double rd = Mathf::Sqrt(Mathf::Pow(r0 + hd, 2) - Mathf::Pow(r0 * Mathf::Cos(input.GetElevation()), 2)) - r0 * Mathf::Sin(input.GetElevation());
            double rw = Mathf::Sqrt(Mathf::Pow(r0 + hw, 2) - Mathf::Pow(r0 * Mathf::Cos(input.GetElevation()), 2)) - r0 * Mathf::Sin(input.GetElevation());

            auto GetAki = [&input, r0](double hi)
            {
                double ai = -Mathf::Sin(input.GetElevation()) / hi;
                double bi = -Mathf::Pow(Mathf::Cos(input.GetElevation()), 2) / (2 * hi * r0);

                std::vector<double> ans(9);
                ans[0] = 1;
                ans[1] = 4 * ai;
                ans[2] = 6 * ai * ai + 4 * bi;
                ans[3] = 4 * ai * (ai * ai + 3 * bi);
                ans[4] = Mathf::Pow(ai, 4) + 12 * ai * ai * bi + 6 * bi * bi;
                ans[5] = 4 * ai * bi * (ai * ai + 3 * bi);
                ans[6] = bi * bi * (6 * ai * ai + 4 * bi);
                ans[7] = 4 * ai * Mathf::Pow(bi, 3);
                ans[8] = Mathf::Pow(bi, 4);

                return ans;
            };

            auto akd = GetAki(hd);
            auto akw = GetAki(hw);

            auto GetDeltai = [](const std::vector<double> &aki, double ri, double Ni)
            {
                double tmp = 0;
                for (int i = 1; i <= 9; i++)
                {
                    tmp += aki[i - 1] / i * Mathf::Pow(ri, i);
                }

                double ans = 10E-6 * Ni * tmp;
                return ans;
            };

            double deltad = GetDeltai(akd, rd, Nd);
            double deltaw = GetDeltai(akw, rw, Nw);

            return deltad + deltaw;
        }

        static DryWet NMF(const TopoCentricPolarCRD &input, double latitude, int curDoy)
        {
            int t0 = 28;
            DryWet ans(6);
            // 计算干分量
            {
                double Pfaita, Pfaitb, Pfaitc;

                if (latitude < 15)
                {
                    Pfaita = Average[a][Angle15] + Amplitude[a][Angle15] * Mathf::Cos(2 * Mathf::Pi * (curDoy - 28) / 365.25);
                    Pfaitb = Average[b][Angle15] + Amplitude[b][Angle15] * Mathf::Cos(2 * Mathf::Pi * (curDoy - 28) / 365.25);
                    Pfaitc = Average[c][Angle15] + Amplitude[c][Angle15] * Mathf::Cos(2 * Mathf::Pi * (curDoy - 28) / 365.25);
                }
                else if (latitude > 75)
                {
                    Pfaita = Average[a][Angle75] + Amplitude[a][Angle75] * Mathf::Cos(2 * Mathf::Pi * (curDoy - 28) / 365.25);
                    Pfaitb = Average[b][Angle75] + Amplitude[b][Angle75] * Mathf::Cos(2 * Mathf::Pi * (curDoy - 28) / 365.25);
                    Pfaitc = Average[c][Angle75] + Amplitude[c][Angle75] * Mathf::Cos(2 * Mathf::Pi * (curDoy - 28) / 365.25);
                }
                else
                {
                    int faii = latitude / 15;

                    Pfaita = Average[a][faii] + (Average[a][faii + 1] - Average[a][faii]) * (latitude - faii * 15) / 15 +
                             (Amplitude[a][faii] + (Amplitude[a][faii + 1] - Amplitude[a][faii]) * (latitude - faii * 15) / 15) * Mathf::Cos(2 * Mathf::Pi * (curDoy - 28) / 365.25);
                    Pfaitb = Average[b][faii] + (Average[b][faii + 1] - Average[b][faii]) * (latitude - faii * 15) / 15 +
                             (Amplitude[b][faii] + (Amplitude[b][faii + 1] - Amplitude[b][faii]) * (latitude - faii * 15) / 15) * Mathf::Cos(2 * Mathf::Pi * (curDoy - 28) / 365.25);
                    Pfaitc = Average[c][faii] + (Average[c][faii + 1] - Average[c][faii]) * (latitude - faii * 15) / 15 +
                             (Amplitude[c][faii] + (Amplitude[c][faii + 1] - Amplitude[c][faii]) * (latitude - faii * 15) / 15) * Mathf::Cos(2 * Mathf::Pi * (curDoy - 28) / 365.25);
                }

                ans[0] = Pfaita;
                ans[1] = Pfaitb;
                ans[2] = Pfaitc;
            }

            // 计算湿分量
            {
                double pfaita, pfaitb, pfaitc;
                if (latitude < 15)
                {
                    pfaita = WetCorr[a][Angle15];
                    pfaitb = WetCorr[b][Angle15];
                    pfaitc = WetCorr[c][Angle15];
                }
                else if (latitude > 75)
                {
                    pfaita = WetCorr[a][Angle75];
                    pfaitb = WetCorr[b][Angle75];
                    pfaitc = WetCorr[c][Angle75];
                }
                else
                {
                    int faii = latitude / 15;
                    pfaita = WetCorr[a][faii] + (WetCorr[a][faii + 1] - WetCorr[a][faii]) * (latitude - faii * 15) / 15;
                    pfaitb = WetCorr[b][faii] + (WetCorr[b][faii + 1] - WetCorr[b][faii]) * (latitude - faii * 15) / 15;
                    pfaitc = WetCorr[c][faii] + (WetCorr[c][faii + 1] - WetCorr[c][faii]) * (latitude - faii * 15) / 15;
                }

                ans[3] = pfaita;
                ans[4] = pfaitb;
                ans[5] = pfaitc;
            }

            auto GetM0iE = [](double E, double ai, double bi, double ci)
            {
                return (1 + (ai / (1 + (bi / (1 + ci))))) / (Mathf::Sin(E) + (ai / (Mathf::Sin(E) + (bi / (Mathf::Sin(E) + ci)))));
            };

            double md = GetM0iE(input.GetElevation(), ans[0], ans[1], ans[2]) + (1 / Mathf::Sin(input.GetElevation()) - GetM0iE(input.GetElevation(), HeightCorrection[0], HeightCorrection[1], HeightCorrection[2])) * input.GetRange() / 1000;

            double mw = GetM0iE(input.GetElevation(), ans[3], ans[4], ans[5]);

            return {md, mw};
        }
    };
}